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ABSTRACT 

The axisymmetric form of the hydrodynamic equations within the smoothed parti- 
cle hydrodynamics (SPH) formaHsm is presented and checked using ideahzed scenarios 
taken from astrophysics (free faU collapse, implosion and further pulsation of a sun-like 
star), gas dynamics (wall heating problem, collision of two streams of gas) and inertial 
confinement fusion (ICF, -ablative implosion of a small capsule-). New material con- 
cerning the standard SPH formalism is given. That includes the numerical handling of 
those mass points which move close to the singularity axis, more accurate expressions 
for the artificial viscosity and the heat conduction term and an easy way to incorpo- 
rate self-gravity in the simulations. The algorithm developed to compute gravity does 
not rely in any sort of grid, leading to a numerical scheme totally compatible with the 
lagrangian nature of the SPH equations. 
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1 INTRODUCTION 

Since the s moothed-particle hydr o dynam ics method was 
proposed bv lGingold fc Monaghanl fToTj) and lLucvl (|l977l ') 
thirty years ago it has become a powerful and very successful 
method, which is routinely used worldwide to model systems 
with complicated geometries. In particular the combination 
of the SPH method with hierarchical algorithms, addressed 
to efficiently organize data and calculate gravity, have led to 
a robust and reliable schemes able to handle complex geome- 
tries in three dimensions. In spite of these achievements little 
effort has been put in the development of axisymmetric ap- 
pUcations of the method. However, a large body of interest- 
ing applications has to do with systems which have axisym- 
metrical geometry such as accretion discs, rotating stars or 
explosive phenomena (i.e. novae or supernovae events) pro- 
vided the ignition takes place in a point-like region at the 
symmetry axis. The interaction between the gas ejected dur- 
ing a supernova explosion and the circumstellar matter has 
been extensivel y simulated us i ng cy lindric coo rdinates either 
with eulerian (jeiondin et all 1 19961 ) or SPH jVelarde et all 
I2OO6I ) codes. In Inertial Confinement Fusion (ICF) studies 
the collapse of the deuterium capsule induced by laser ab- 
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lation can be also approximated using two-dimensional hy- 
drodynamics. Axisymmetrical hydrocodes are also useful to 
conduct resolution studies of three-dimensional codes just by 
imposing an initial configuration with the adequate symme- 
try. In these cases a comparison between the results obtained 
using the 3D-hydrocode and the equivalent, although usually 
better resolved, 2D simulation could help in numerical con- 
vergence studies. Several formulations of different complex- 
ity have been proposed to handle with the axisymmetric ver- 
sion of SP H. While the fir s t algo rithms were in general very 
crude, i.e. iHerant fc Ben3 (| 19921 ) and references therein, ne- 
glecting the so called hoop-stress te rms, there were also 
more consistent SPH approximations, (|Petscheckfc Liberskvl 
:199 sj)- Recently more co mplete formulations were given by 
Bro okshawl (|2003l ) and bv lOmang et all l|2006l ). In particular 
the scheme devised by Omang et al. includes not only the 
hoop-stress terms but also a consistent treatment of the sin- 
gularity line defined by the symmetry axis. The algorithm 
proposed by these authors relies in the use of interpolative 
kernels especially adapted to the cylindrical geometry. Al- 
though the resulting scheme is robust, being able to success- 
fully pass several test cases, it has the numerical inconve- 
nient that an elliptical integral has to be solved numerically 
at each integration step for each particle. 

In this paper we propose a new formulation of the 
axisymmetric SPH t echnique which p reser ves many of 
the vi rtues stated by iBrookshawl (|2003l ) and lOmang et al.l 
l|2006l ). being also able to handle the singularity axis in an 
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efficient way. In our scheme there is no necessity to modify 
the interpolating kernel (we use the cubic spline polyno- 
mial). Instead, taking advantage of the symmetry proper- 
ties across the Z-axis, we found several analytical correction 
factors to the physical magnitudes. These correction factors 
only affect particles moving at a distance from the Z-axis 
lesser than 2h, being h the smoothing-length parameter. 
For the remaining particles the SPH equations are identi- 
cal to those given in lBrookshawl (2003 ) . In addition, we also 
give expressions for the artificial viscosity and a new heat 
conduction algorithm, which takes into account the hoop- 
stress contribution. Finally, we propose an original method 
to handle gravity within the SPH paradigm which is consis- 
tent with the gridless nature of that particle method. The 
text is organized as follows: the basic fluid equations written 
in axisymmetric coordinates, and corrected from axis eflects 
when necessary, are provided in Section 2. In Section 3 we 
add useful physics to these equations consisting of an ar- 
tificial viscosity term to handle shock waves (Section 3.1), 
an expression for thermal conduction (Section 3.2) and an 
algorithm addressed to calculate self-gravity within the ax- 
isymmetric hypothesis (Section 3.3). Section 4 is devoted to 
describe and discuss five test cases aimed at validating the 
proposed scheme. Finally the main conclusions of our work 
as well as some comments about the limitations of the de- 
veloped scheme and prospects for the future are presented 
in Section 5. 



2 FLUID EQUATIONS IN THE 
AXISYMMETRIC APPROACH. 

An elegant approa ch to the ax i symm etric Euler fiuid equa- 
tions was given bv lBrooksha^ (|2003l ) who derived the SPH 
form of these basic equations using the minimal action prin- 
ciple (see Monaghan 2005 and references therein for the his- 
tory of variational SPH). The resulting expressions for mass, 
momentum and energy conservation written in the cylindri- 
cal coordinate system s — (r, z) are: 

N 

r]i = ^rnjWij (1) 

', = 2.f-2.£m,(iV. + £i^)^ 

- ^7 ' ' (2) 

^^_2^:^+2^:^ym,(vi-Vj).D. (3) 

where rji — 2-nripi is the two-dimensional density of the i- 
particle, v — {vr,Vz) its velocity, Wij = Wij{ ^^\^''^ ) the 
interpolating kernel and the remaining symbols have their 
usual meaning. The differential operator D = {-§p, -§^) is 
a cartesian operator written in the (r, z) plane of cylin- 
dric coordinates. The first terms on the right in the r- 
component of equation (2) and equation (3) are called the 
hoop-stress terms. They are geometrical terms which arise 



from the nabla operator written in cylindrical coordinates. 
The smoothing-length parameter, h, is usually taken as the 
local resolution of the SPH. Equations (1), (2) and (3) 
along with the adequate equation of state (EOS) are the 
starting point to carry out numerical simulations of fiuid 
systems with axisymmetric geometry. As boundary condi- 
tions we consider refiective ghost particles across the Z-axis 
and open flow at the outer limits of the system. For the 
i-particle with coordinates {ri,Zi) and velocity {ri,Zi) it is 
defined its reflective fc-particle by taking (rfc, Zk) = {—f'i, Zi), 
{vk, Zk) = {—ri, Zi) and rak ~ rrii. Thus position, velocity as 
well as other magnitudes such as density or internal energy 
of refiective particles are updated at each step not using 
the SPH equations but from the evolution of real particles. 
Even though the use of reflective boundary conditions is not 
strictly necessary in axisymmetric geometry it is useful to 
correctly represent the density and its derivatives near to 
the singularity axis. 

2.1 Correction terms close to the singularity axis 

One of the major causes of inaccuracy of axisymmetric SPH 
is the treatment of particles that get close to the symmetry 
axis. Unlike spherically symmetric systems where there is 
only one singular point, just at the center, here there is a 
singular line at r = 0. A good treatment of particles mov- 
ing close to the Z-axis is especially relevant for implosions 
such as the collapse of a self-gravitating system or in in- 
ertial confine ment fusion studies. As far as we know only 
lOmang et al. I 12006,) 

have consistently addressed this prob- 
lem. In that paper the interpolation kernel is modified ac- 
cording to the particular geometry of the system, spherical 
or cylindrical. The resulting scheme does not have singular- 
ity problems when the particles approach the axis. However 
the resulting kernel does not have an analytical expression 
and must be calculated at any step for each particle using a 
numerical in tegration. Rec e ntly, u seful fitting formulae were 
proposed bv lOmang et al.| |20o3) but still involving a large 
number of operations which slows the calculation. 

An alternative way to handle with the symmetry axis, 
without modifying the basic SPH scheme given above, is 
to calculate correction terms to equations (1), (2) and (3), 
which become significant only close to the Z-axis. These 
correction factors arise because of the limited capability of 
standard kernels to interpolate accurately non linear func- 
tions. In the particular case of the two-dimensional density, 
rj = 27r|r|p, its profile is not longer linear in the axis neigh- 
bourhoods due to the presence of refiective particles. Usually 
the errors are small and the interpolation is precise to second 
order in h. Unfortunately, close to the symmetry axis errors 
grow and density and other physical magnitudes are not well 
reproduced, as it will be shown below. The detailed deriva- 
tion of these corrections are given in Appendix A, leading 
to the following equations: 

N 

=ym,miX/{ = r,. x/1 (4) 

where rji is the new, improved, two-dimensional density and 
fl is a correction factor which, for the cubic-spline kernel, 
reads (Appendix A): 



Axisymmetric smoothed particle hydrodynamics with self-gravity 3 



2.5 



1.5 - 



o 
■t-i 
o 
(0 

c 
o 

o 

0) 

1 

o 
o 



X 

< 



0.5 



-1 — I — I — r- 



n — I — I — T" 



X 



f, 



0.5 



1.5 



Figure 1. Plot of the correction factors /i and /2 and their 
derivatives as a function of C = r/h. The cubic spline interpolator 
kernel was assumed to compute /i and /2 (Appendix A). 
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being (i = Ti/hi. A plot of /i(C) and its first derivative 
is given in Fig. 1. Hereafter a hat will be placed over any 
corrected, and therefore, "true" magnitude. 

At this point it is useful to introduce a couple of alge- 
braic rules: 

A, 



(6) 



(7) 



which are valid as long eis, close to the axis, the magnitude 
A has a weak dependence in the r-coordinate. In particular, 
setting ^ = 1 in equation (7) leads to equation (4). 

To find the influence of the above correction factor in 
the momentum equation, equation (2), its is better to write 
the components of that equation in differential form. For the 
radial component we have: 



1 dP 
p dr 



2nr dP 
rf dr 



P 2nrP 8^ 
dr 



V 



d f 2-KrP\ 
dr \ rj J 



When r — » the particle approaches the Z-axis where sym- 
metry enforces the last term on the right to vanish for spher- 
ically symmetric kernels. The first term on the right is the 
hoop-stress term which should be exactly balanced by the 
central term when r — > 0, because the acceleration should be 



zero at the symmetry axis. Nevertheless, this does not hap- 
pen unless the correction factor /i(C) is taken into account 
during the calculation of the gradient of rj: 



dy _ d{MOv) _ dr, dMC) 

dr dr ■'^^^'dr^ dr 



(9) 



Similarly, for the Z-component in the momentum equa- 
tion: 

.. _ IdP _ 2irr dP _ 2nrP drj d f 27rrP \ 
p dz rj dz rp dz dz \ rj J 

where now we simply have: 



d-z~^'^^^d-z 



(11) 



The resulting momentum components in discrete SPH 
form can be obtained from equations (6) to (11): 



27r5- - 2n 



N 

E 



( 2nPrA dfiiCi) 



yVi X fi J dr. 



(12) 



-2nJ2^ 

3 = 1 



dzi 



(13) 



The value of pressure in equations (12) and (13) must be 
computed through the EOS using the 'corrected' 3D-density, 

p = ?j/{2T:r). 

Note that close to the axis the summatories in equa- 
tions (12) and (13) are non symmetric with respect any pair 
of particles. Still the r-component of the momentum is con- 
served, even for those particles with r < 2h, because of the 
imposed reflective boundary conditions. In the Z-axis the 
momentum is not exactly conserved within the band de- 
fined by r < 2h. In general such loss is small, affecting a 
tiny amount of mass. In most of the tests described below 
total momentum was very well preserved. 

A similar approach can be worked out to improve the 
energy equation: 



du 
It 

where: 



P 

-2n— Vr 



+ 2ix 



Pr drj 
'Wdi 



(14) 



(8) 



drj 
lit 



drj d{riVr) 
dr dr 



+ 



drj 8{r,Vz 
' dz dz 



(15) 



being rj=r]X fi; (rjVz) = {r]Vz) x /i and {r]Vr) = {rjVr) x /2. 
The factor is (Appendix A): 



4 D. Garcia-Senz et al. 



f [f|C^ + |0-liC?+e|jC*]-' 0<0<1 

r L/'"^ I 16 A-l „ 1 I 8 A _ 1^2 I _La3_ 

L 45 15^» 3 ' 15 

tmCA for K < 2 



(16) 



1 > 2 

The resulting energy equation in discrete SPH form is: 



dui 



at rii r/i dt 



(17) 



where: 



J = l 

JV 



9ri 



i=i 



i=i 



(18) 



Note that fl, f2 and their derivatives are only func- 
tion of the current radial coordinate of the particle and its 
smoothing length. Thus, they can easily be computed and 
stored in a vector without introducing any significant com- 
putational overload. Examples about the use of the axisym- 
metric fluid equations with axis corrections, equations (4), 
(12), (13) and (17) in practical situations will be provided 
and discussed in Section 4. 



3 ADDING PHYSICS: SHOCKS, THERMAL 
CONDUCTION AND GRAVITY 

3.1 Artificial viscosity 

As in three dimensions the treatment of shocks in the ax- 
isymmetric approach also relies in the artificial viscosity for- 
malism. Nonetheless, there are a variety of artiflcial viscosity 
algorithms suited f or SPIf to choose at. We have adapted the 
standard recipe by iMonaghan fc Gingoldl (|l983l ) to the pe- 
culiarities of the 2D-axisymmetric hydrodynamics. In that 
approach the artificial viscosity gives rise to a viscous bulk 
and shear pressure which is added to the normal gas pres- 
sure only in those regions where the particles collide (in the 
SPH sens43)- In three dimensions it is defined a magnitude, 

n„: 



'-'-1. 



Vij • rij < 



(19) 







which is closely related to the viscous pressure. In equation 
(19) ao and /3o are constants of the order of unity, pij and 



Cij are the average of density and sound speed of i and j- 
particles, and fiij is: 



hi 



r + h"' 



(20) 



here Vij — (v^ — Vj), Vij = (r^ — r^) and v = Q.lh avoids 
divergences when rij —> 0. The scalar magnitude Ilij has 
a linear and a quadratic term. The linear component mim- 
ics bulk and shear viscosity of fluids whereas the quadratic 
one is important to avoid particle interpenetration in strong 
shocks. 

The easiest way to write the contribution of the artifl- 
cial viscosity to the momentum and energy equations in an 
axisymmetric code is by changing the mass of the particles 
according to their distance to the Z-axis: m — > if^. The 
viscous acceleration becomes: 



3D 



N 

i=i 



(21) 



Taking pij — '7ij/(2'rrij) in equation (19) the explicit 
dependence on rij in the viscous acceleration formula is 
removed. 



j=i 



(22) 



where li^^J is 



n 



(1) 



Vij ■ Sij < 







being pij-. 



Sij + V' 



(23) 



(24) 



Expressions (22), (23) and (24) account for the cartesian 
part of viscosity. It has been shown (i.e. Monaghan 2005) 
that, in the continuous limit, these expressions become the 
Navier-Stokes equations provided shear and bulk viscosity 
coefficients are taken rj^is = paohc/8 and (^vis ~ 5r;i,is/3 re- 
spectively. However, in cylindric geometry the stress tensor 
which appear in the Navier-Stokes equations also includes 
a term proportional to velocity divergence through the so 
called second viscosity coefficient. Therefore an extra term 
containing the magnitude (vr/r) must be added to H'j'' to 
account for the convergence of the fiux towards the axis. 

(2) 

This new term, H^^. , should fulfill a few basic requirements: 
1) Far enough from the axis it should be negligible, 2) must 
vanish for those particles with v,. > 0, 3) for homologous 
contractions the viscous acceleration related to that term 
should be negligible, 4) it should be symmetric with respect 
particles i and j to conserve momentum. An expression sat- 
isfying the four items is: 







u^ .and Vr <Q 



other cases 



(25) 



^ In the SPH method the so called particles could be looked as 
finite spheres of radii 2h 



where qij — 0.5 ( ' ''' -| — - — -). The resulting viscous accel- 
eration is: 
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(26) 



where 



n 



(1) 



(2) 

Note that for a homologous con- 



traction Vr cc r meaning q ~ const., hence the gradient of 
<^?7n'^^^ vanishes fulfilling the third requirement above. In 

(2) 

diverging shocks 11^^ vanishes and only the cartesian part 
of viscosity matters. Therefore constants qo,/3o should re- 
main close to their standard values qq — 1 and /3o — 2. All 
simulations presented in this paper were carried out using 
ai = Qo = 1, I3i — /3o = 2. In converging shocks, however, 
the effect of 11^^ is to increase artificial viscosity introducing 
more damping in the system. 

Equation (26), has the advantage that is formally simi- 
lar to that used in three-dimensions. Therefore one can ben- 
efit from the well known features of the artificial viscosity 
in 3D, which can be directly translated to the axisymmet- 
ric version (several useful v ariations of equation (19) can 
be found in iMonaghan 20051 ). In particular, it is straightfor- 
ward to write the corresponding energy equation: 



' dui 



1 ^ 

-^mjn-/'(vi 



Vj)-D,m, 



(27) 



which has to be added to the right hand of equation (3). 



3.2 An approach to the conduction term. 

The differential equation describing the evolution of the 
specific internal energy due to conductive or diffusive heat 
transfer is: 



'du\ 1„ , „ , 

— ) =-v-(Kvr) 

^ ox J cond P 



(28) 



being k the conductivity coefficient which in turn is a func- 
tion of the local thermodynamic state of the material. The 
main difficulty to write equation (28) in a discrete equation 
suitable to SPH calculations is the existence of a second 
derivative. It is well known that second and higher order 
derivatives often pours a lot of numerical noise in disordered 
systems. A way to avoid that shortcoming is to reduce one 
degree the order of the der ivative by t a king t he integral ex- 
pression of equation (28 ). iBrookshaw] 1 1985h . It has been 
shown (see, for example. lJubelgas et al. 1 2004 )) that in 3D- 
cartesian coordinates the following expression allows to ap- 
proximate a second derivative using only the first derivative 
of the interpolation kernel: 



N 



P3 



r^J ■ ViWij 



(29) 



where Y represents any scalar magnitude and r^j = Vi — Vj . 
An useful algebraic relationship which allows to write the 
heat transfer equation in SPH notation is: 



V • (kVT) = - [V^(kT) - TV^K + kV^t] 



(30) 



Combining equations (29) and (30) the usual form of the 
heat transfer equation used in 3D-SPH studies is obtained. 
A mathematical expression for the heat transfer e quation in 
axisymmetric SPH was given bv lBrooksha^l|2003l ). However 
in the derivation of equation (29) there was not taken into 



account the ~ 1/r term which naturally arises whenever the 
divergence of a vector is estimated in cylindric coordinates: 



V • (VF) 



r dr \ dr J dz \ dz J 



(31) 



It will be shown in Section 4.1 how the inclusion of the 
first term on the right of equation (31) improves the quality 
of the results in a particular test. 

To obtain the adequate numerical approximation to 
equation (28) in axisymmetric SPH we first make use of 
equation (30): 



2p ^ ' 



(32) 



Then, using equation (31) to develop each term inside the 
brackets and using the 2D-operator D = {■§^, -§^) instead 
of V the following expression is obtained: 



^dui\ 

\ dt / cond 



1 j 1 djKTj) Tj dm ^ K,i ar, 

2p 1 Ti dn Ti dri Ti dvi 



+ 



N 

1 



Pi 



^^i ~\~ f^j 
} ~ 



(33) 



The magnitudes inside the brackets are the new terms 
related to the hoop-stress. In order to calculate the deriva- 
tives we make use of one of the golden rules of SPH 
( Monaghan 2005 ): 



da 



dr. 



Vj 



dri 



(34) 



After a little algebra the following expression is obtained: 



\ at J cond ^ 



mm 



dn 



i=i 



m-T,)4-D,VK,, (35) 



The presence of (Ti — Tj) in the equation ensures that 
there is not heat flux between different parts of an isother- 
mal system. Note that the presence of the multiplier in 
the second term on the right side of equation (35) does not 
ensure complete conservation of the heat flux. However the 
total energy losses in the numerical test below simulating 
a thermal wave evolution were negligible. Of course equa- 
tion (35) can be symmetrized by taking the arithmetical 
mean fij instead of but in that case the evolution of the 
thermal wave was not so well reproduced. On another note 
equation (35) is compatible to the second principle of ther- 
modynamics in the sense that heat always flows from high to 
low temperature particles. To demonstrate this let us take 
a pair of particles i and j so that Ti > Tj. The second term 
on the right becomes negative because the scalar product 
Sij • HiWij is always negative. The first term on the right is 
negative for < rj and positive for > rj . As the heat flow 
from particle i must be negative the only trouble could come 
if ri > rj. Nevertheless, even in that case heat flux is still 
negative if rijiji — rj) > 1/2 because the sign of the second 



term in equation (35) prevails. Thus ri/{ri 



> 1/2 is 
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Figure 2. Sketch of the coordinate system and notation used to 
describe gravity. 



a sufficient condition the get the right sign of heat flux be- 
tween any pair of particles. As r; — r j ~ /i such condition is 
fuUfiled in a large domain of the system. The exception could 
be the axis vicinity where — > 0. Nevertheless in that case 
symmetry enforces the heat flux to be negligible. Therefore 
we expect a good behaviour of the heat flux arrow although 
there are not excluded marginal violations if resolution is 
poor and strong heat fluxes were present close to r = 0. A 
way to ensure complete compatibility to the second princi- 
ple of thermodynamics is to make zero the flux between any 
pair of particles violating such principle. 

3.3 Self-Gravity 

Current 2D-hydrocodes often handle gravity by solving the 
Poisson equation or, if the system remains nearly spheri- 
cal, by simply computing the enclosed lagrangian mass in a 
sphere below the point and using the Gauss law. Methods 
based on the Poisson solvers have proven very useful to find 
gravity in eulerian hydrodynamics where the same grid used 
to compute the motion of the fluid elements can be used in 
the calculation of gravity. However they have the difficulty 
to set suitable outer boundary conditions owing to the infi- 
nite range of the gravitational force. In lagrangian gridless 
methods such as SPH it is better to use the direct interaction 
among mass particles themselves to calculate gravity. The 
evaluation of the gravitational force through direct particle- 
particle interaction leads to an A'^'^ scheme that makes the 
computation feasible only for a limited number of particles. 
When N is high, as it is frequent in three-dimensional cal- 
culations, one has to rely in approximat e schemes such as 
those based in hierarchical-tree methods. iHernguist fc Katj 
However hierarchical methods do not work efficiently 
in the 2D-axisymmetric approach because what we call par- 
ticles are in fact rings of different size. Usually the ratio 
between the radius of these rings and the distance to the 
point where the force needs to be computed is too large to 
permit the multipolar approach to evaluate the gravitational 
force. Fortunately, the good resolution usually achieved in 
2D using a moderate number of particles makes the direct 
calculation affordable. 

According to Fig. 2 the gravitational force per unit of 



mass in a point P of coordinates (0,j/p,Zp) (being Z the 
symmetry axis) due to the ring is: 



= G 



pRdtp 



(7?2 + s'2_22/pi?sin(^)3/2 

[(i?sin - Vp) i + {z - Zp)\i\ (36) 



where the meaning of the symbols is that shown in Fig. 2. 
The gravitational force acting onto the i-particle can be eas- 
ily written as: 



gi = 



^ 2^ [R] -f )3/2 



[(i?,/i -y,/2) j + (2j -2«)/2) k] 



(37) 



where a'^ = + {zj — Zi)'^ , and rrij = 2TvRjPj is the mass of 
the particle associated to the j-ring. Integrals Ii and I2 are 
defined as follows: 



sin 93 d(f 



(38) 



(39) 



(40) 



(1 — T sin 1^)3/2 
d(p 

(1 — T sin iy5)3/2 
where the parameter t is: 

^ 2y,Rj 
R]+sf 

Although the elliptical integrals /i and I2 can not be solved 
analytically they can be tabulated as a function of the pa- 
rameter T given by equation (40). It is straightforward to 
show that the value of r is always inside the interval Te[0, 1), 
although for r ^ 1 the integrals I\ and I2 become divergent. 
Therefore the gravity force can be calculated in an efficient 
way using the following recipe: 

1) Build a table for I\ and I2 as a function of r. A table 
with 10^ rows with values ^ r ^ 0.9999 evenly spaced is 
sufficient. 

2) To increase the speed do not interpolate from that 
table but take just the row which is closest to the actual 
value of r calculated using equation (40). 

3) Note that parameter r is symmetric with respect any 
pair of particles, Ti = Tj, thus /i(ri) = Ii{Tj) and hin) = 
72 (tj). Therefore only half of the interactions have to be 
calculated. 

If the algorithm is optimized the scheme is able to pro- 
vide the exact value of the gravity for several dozens of thou- 
sand particles. In many applications using ~ 5 10* particles 
in two-dimensions is enough to guarantee a good resolution. 

Another physical magnitude of interest is the gravita- 
tional potential at the position of the i-particle. It is easy to 
show that the contribution of the j- ring to the gravitational 
potential is: 



V, = 



G 



2n (i?^ + sf)i/2-^3 



where, 



/3 = 



d(p 



(J {1 — T sin ip)^^^ 



(41) 



(42) 
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Again the same recipe given above to compute the force 
can be used to efficiently calculate Vi. On the other hand 
the exact computation of the gravitational potential allows 
to calculate the gravitational force by taking the gradient of 
V at any point. 



gi = -i'DV)i 



(43) 



A more suitable form for SPH calculations can be ob- 
tained using: 



(44) 



which, according to equations (7) and (9), leads to the fol- 
lowing discrete equation: 



{-DV)i = iij2mj{Vi- Vj)DiWij + 



m 



Vi 



^mj{Vi - Vj)Wij Ur 



(45) 



where /i is the corrective term given by equation (5) and 
Ur the unit vector in the r-axis. 

This second route to calculate the gravitational force 
is computationally more efficient than evaluating equation 
(37) because the gradient of the potential is a local quantity 
which can be calculated in the same part of the algorithm 
devised to compute the density or other magnitudes in the 
hydrocode. It has the additional advantage that the result- 
ing force is smoothed by the SPH interpolation procedure 
avoiding divergences when a pair of particles become too 
close. In Fig. 5 (bottom-right) there are shown the gravity 
profile calculated using equation (45) (filled triangles) and 
the pressure gradient term (continuum line) along a sun-like 
polytropic structure. As it can seen the fit is satisfactory 
except at the surface where the pressure gradient is overes- 
timated. Although using the potential to calculate gravity 
is not as exact as the direct force calculation it is a fac- 
tor two faster because there is a lesser amount of numerical 
operations in the double loop of the gravity routine. 

Needless to say, the simplicity of the proposed scheme 
makes the parallelization of the gravity computational mod- 
ule straightforward. In this case calculations with 10® — 
10^ particles could become feasible even for desktop com- 
puters with multiple core processors. 



3.3.1 Free-fall collapse of homogeneous gas structures. 
Rotation. 

As an initial check of the gravity algorithm resulting from 
equation (45) we have simulated the free-fall collapse of a 
uniform density sphere of mass Mo and radius Rq. It is a 
standard test which has the following analytical solution: 



iff 



where ro is the initial position of the fiuid element and ^ 
t ^ tff. The free-fall time iff is: 



0.8 



0.6 



tt! 0.4 



0.2 



n — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — T" 



^ v,^ Homogeneous rotating 
> cylinder 



.. Homogeneous cyiinder 



Homogeneous sptiere 



t/tff 



Figure 3. Free-fall test. Triangles represent the analytical so- 
lution in each case. The initial radius of the particle was ro = 
2/3 Rq, 20 = in all cases. 



tff = 



2 \2GMo 



(47) 



We built an uniform sphere with Mo=l Mq and 
Ro=R0 filled with N = 15396 particles settled in a square 
lattice. Gas pressure and artificial viscosity were set to zero 
so that the structure collapsed under gravity force. After- 
wards the implosion was followed until the elapsed time was 
close to tff. Although restricted to spherical symmetry the 
free-fall test is quite demanding because the evolution is 
highly non linear, allowing for a good check of both the 
gravity module and the integration scheme (see next section 
for details). In Fig. 3 there is represented the evolution of 
a particle initially located at ro = 2/3 Rq. As we can see 
its evolution is in good agreement to the analytical solution 
given by equation (46). 

A question of great interest in astrophysics is the ca- 
pability of axisymmetric SPH codes to handle rotation. 
The topic is far from trivial because in general it involves 
transport of angular momentum via viscous coupling. Even 
though a complete answer to that question is beyond the 
scope of the present work there is a particular case that can 
be handled with the present scheme: the fast implosion (or 
expansion) of a self-gravitating rotating cloud. If the char- 
acteristic dynamical time is much shorter that the viscous 
coupling time we can impose angular momentum conserva- 
tion around the symmetry axis to solve this problem. The 
strategy is simply to add the centrifugal force which arises 
from finite angular momentum to the r-componcnt of grav- 
ity. As an example we have simulated the collapse of a slen- 
der rotating cylinder of gas with uniform density centered at 
the coordinate origin. The initial conditions are specified by 
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the mass M^yi the radius Rcyi and the length of the cyhn- 
der Zcyi- If we suppose rigid rotation the specific angular 
momentum of a mass point is given by Lz(ro) = wo Vq be- 
ing Wo the initial angular velocity of the cylinder and ro the 
position of the fluid element at the initial time to- Angular 
momentum conservation demands L^{r,t) — Lzira) so that 
the centrifugal force is Fc = Ll{ri)) / r'^ (t) which is added 
to the first component of gravity, calculated using equation 
(45), to obtain an effective gravity value. As in the case 
of the spherical collapse we have taken M^yi = Mq and 
Rcyi = Rq while the cylinder length was Zcyi = 20Rcyi. 
The evolution of mass points were followed assuming two 
values for the angular momentum: 1) zero angular momen- 
tum and 2) Lz{ro) = (0.5gor^) where ro = 2/3 Rcyi for 
which centrifugal force was a half of gravitational force at 
that position. In this case we have used a larger number of 
particles, A'' = 70000 to reduce boundary effects. 

On the other hand, assuming Zcyi » Rcyi and us- 
ing the Gauss-law it is possible to work out an analytical 
approach to the above scenario. The acceleration equation 
becomes: 



r = -sro— + -5- 



(48) 



where Lo = iz(ro) and go is the value of gravity at r = ro. 
According to the Gauss-law go — 2GMcyi ro/{Rcyi Zcyi), 
but it is better to take go directly from the SPH simulation 
to ensure identical initial conditions in both calculations. 
The solution of equation (48) is: 



t{r) 



dr' 



-2goro In 



- + ^0 f ^ 



(49) 



which has to be solved numerically once ro and Lo are spec- 
ified. 

Fig. 3 depicts the evolution of a particle initially located 
at ro = 2/3 Rcyi,Zo — without and with initial angular 
momentum as well as that obtained using the the analytical 
approach, equation (49). The case with zero angular mo- 
mentum led to the free fall of the mass element which was 
less violent than the spherical case owing to the lower initial 
density in the cylinder. As it can be seen in Fig. 3 the agree- 
ment between analytical and SPH results is not as good as 
in the spherical case for t > tff. This is not surprising be- 
cause boundary effects at cylinder edges progressively affects 
gravity at current particle test position and its evolution is 
very sensible to small variations of gravity force. 

When angular momentum was added to the cylinder at 
to = s the implosion of the structure slowed down. As com- 
mented above the amount of angular momentum was chosen 
to get a centrifugal force contribution at ro = 2/3 Rcyi equal 
to a half of gravity force at that position. As we can see in 
Fig. 3 the result of the simulation is in better agreement 
to the analytical approach than the pure free-fall case. This 
is because the evolution is driven not only by gravity and 
errors due to the finite size of the rotating cylinder are not 
affecting so much the outcome as in the non rotating case. 
At t ~ 5 t// the fluid element begins to be centrifugally 
sustained, in good match to the analytical estimation. 



4 TEST CASES 

In this section we describe and discuss in detail five test cases 
addressed to validate the computational scheme developed 
in the precedent sections. The first test involves the prop- 
agation of a thermal discontinuity born at the symmetry 
axis. It is a well known calculation which has an analytical 
solution to compare with. Three calculations: the gravita- 
tional collapse of a polytrope, the implosion of a homoge- 
neous capsule, and the wall shock problem deal with implo- 
sions of spherically symmetric systems. Although, at first 
glance, such constraint looks too restrictive in fact it is not, 
because the spherical symmetry is not a natural geometry 
for axisymmetrical systems described with cylindrical co- 
ordinates. Any deviation from the pure spherical symmetry 
during the implosion will trigger the growth of hydrodynam- 
ical instabilities. Therefore the preservation of the symmetry 
during the calculation is an important feature to be added to 
mass, momentum and energy conservation. Finally the last 
test was devoted to simulate the collision of two bubbles of 
fluid along the Z-axis. 

The initial models were calculated by mapping the 
spherically symmetric profiles into a 2D distribution of parti- 
cles settled in a square lattice. The mass of the particles was 
conveniently adjusted to reproduce the density profile of the 
one-dimensional model. For example, masses proportional to 
the r— coordinate of the particle were taken to obtain models 
with constant density. Reflective boundary conditions using 
specular ghosts particles were imposed across the Z-axis. 
The EOS was that of an ideal gas with 7 = 5/3 and radia- 
tion in the test devoted to the collapse of a polytrope, and 
only gas ideal for the other tests. We have used and adap- 
tive smoothing length parameter h{s, t) which was updated 
at each time step to keep a constant number of neighbours 
n„ = 36 within a circle of radius 2h. The interpolating ker- 
nel was the cubic polynomial spline. Numerical integration 
of SPH fluid equations was performed using a two-step cen- 
tered scheme with second order accuracy. It is worth to men- 
tion that in the free-fall test of a homogeneous sphere the 
velocity was updated using the XSPH variant of iMonaghanI 
(| 19891 ) to avoid particle penetration through the symmetry 
axis. In this particular test mass points at low radius are 
proner to cross Z-axis because pressure and viscous forces 
were artificially set to zero. Therefore any residual value of 
gravity can impel particles located close to the center to- 
wards unphysical r < values. We have updated only the 
r-component of velocity using the following expression: 



n = r, 4- e 27r r^ 5^ (f, - r^) W„' 



(50) 



where r^ is the new, smoother, velocity. The parameter e was 
taken variable: 



eo 







i_(ir 



r ^ 3/1 



(51) 



r > 3/1 



where eo = 0.5. Using XSPH variant to move particles en- 
forces f ^ when r 0, making it difficult for particles to 
cross the axis. 
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4.1 Evolution of a thermal discontinuity 

Lets consider the problem relative to the propagation of a 
thermal wave moving through an uniform static medium. 
This is a well known test, addressed to check the capability 
of the numerical scheme to handle thermal discontinuities. 
The initial model was a sample of 57908 particles evenly 
distributed in a square lattice so that the density was p = 
1 g.cm~^. A (5-like jump in internal energy originates a a 
thermal wave front which evolves according to: 

where c„ is the specific heat and k is the thermal con- 
ductivity. The following set of values were taken: A = 
10^ ergs.cm^/g, uo = 10^ erg.g"^ and c„A; = 1 cm^s"^ The 
initial internal energy profile was that given by equation (52) 
for t = 1 s, which was taken as the initial time (t=0 s) for 
the SPH simulation. The evolution of the thermal signal is 
then controlled by the heat conduction equation, equation 
(35). 

In Fig. 4 (left) there is represented the thermal profile 
at different times. As we can see the coincidence with the 
analytical solution given by equation (52) is excellent. As 
time goes on the peak of the signal and its slope decreases 
due to heat diffusion. The initial discontinuity is rapidly 
smeared out by thermal diffusion and soon a thermal wave 
is born which travels to the right, equalizing the internal 
energy of the system. At t=5 s the profile of the internal 
energy of the gas is already quite flat and the system is not 
far from thermal equilibrium. At that moment total energy 
was conserved up to 5 10~^. In Fig. 4 (right) it is shown 
the evolution of the thermal profile which results when the 
first term on the right side in equation (35) is removed. As 
we can see the evolution is no longer reproduced by the 
SPH calculation. Therefore it is of utmost importance to 
include that term, especially in those calculations dealing 
with strong thermal gradients close to the symmetry axis. 

4.2 Gravitational collapse of a polytrope 

The second test involves a catastrophic, albeit highly un- 

probable, astrophysical scenario. A spherically symmetric 
sun-like polytrope was suddenly unstabilized by removing 
the 20% of its internal energy so that the structure collapsed 
under the gravity force. At some point the collapse is halted 
and an accretion shock forms which manages to eject part of 
the mass of the polytrope. Several episodes of recontraction 
followed by mass loss ensued until the star sets in a new equi- 
librium state. Even though that particular scenario is not 
realistic it contains several pieces of physics of great inter- 
est because accretion shocks and pulsational instabilities are 
ubiquitous in astrophysics. As the initial model has spherical 
symmetry we expect it to be preserved during the implosion 
and further bounce. The conservation of the symmetry is 
a demanding test for multidimensional hydrocodes in those 
cases where there are episodes of strong decelerations. In the 
particular case of axisymmetric hydrodynamics the higher 
numerical noise close to the symmetry ax;is may trigger the 
growth of convective instabilities. An additional advantage 
of considering an spherically symmetric initial model is that 



the evolution calculated with the SPH code can be checked 
using standard lagrangian hydrodynamics in one dimension. 

The initial model was a 1 Mq spherically symmetric 
polytrope of index n = 3. The radius was set equal to 
1 Rq so that the central density was pc — 77 g.cm"'^ . Once 
the ID equilibrium model was built it was mapped to a 2D 
distribution of 51408 particles located in a rectangular lat- 
tice. The mass of the particles were conveniently adjusted to 
reproduce the density profile of the polytrope. In Fig. 5 there 
arc shown the profiles of density and gradient of pressure at 
t=0 s calculated using the 2D-SPH. As we can see the inclu- 
sion of the corrective term /i given by equation (5) in the 
momentum equation is crucial to get good enough profiles 
of these quantities to guarantee the stability of the initial 
model. The polytrope was perturbed by reducing the tem- 
perature everywhere in a 20% of its equilibrium value. Af- 
terwards the evolution was followed with the 2D-SPH from 
the implosion until the first pulsation and compared to that 
obtained by a ID lagragian hydrocode. The main features of 
the model and a summary of the results are shown in Table 
1. 

Soon after the model was destabilized the polytrope 
started to collapse. At t= 960 s a maximum of the cen- 
tral density pmax ~ 192.8 g.cm^'^was reached. A similar 
maximum of pmax = 187.1 g.cm^'^was obtained using one- 
dimensional hydrodynamics. The profiles of density and rar 
dial velocity at different times are depicted in Fig. 6. As we 
can see the evolution calculated in one and two dimensions is 
very similar and, in general, both profiles are in good agree- 
ment. At our last calculated time, t= 1545 s, the shock is 
already brealcing out the surface of the polytrope. Shortly 
after that time some mass is ejected from the surface and, 
as the ID calculation shows, the star embarks in a long pul- 
sational stage until a new equilibrium state is achieved. 

Therefore the numerical scheme was able to handle with 
this scenario. The algorithm devised to calculated the grav- 
ity using equation (45), which relics in the direct interaction 
between rings (Fig. 2), did a good job. The artificial vis- 
cosity module was also able to keep track with the shocks 
although, at some stages, the post-shock region showed a 
small amount of spurious oscillations. These numerical os- 
cillations in the velocity profile are clearly visible in Fig. 6 
at t= 1545 s. In the last three columns of Table 1 there 
is information concerning momentum and energy conserva- 
tion at t= 1545 s, our last calculated model. There was an 
excellent momentum conservation, close to machine preci- 
sion, whereas the conservation of energy was more modest, 
~ 0.5%. On the other hand the spherical symmetry was also 
quite well preserved during the calculation. On the negative 
side we find that the two-dimensional calculation was sys- 
tematically delayed with respect its one-dimensional coun- 
terpart. The relative shift in time remained approximately 
constant, around 1, 5%, during the evolution. For the sake 
of clarity the elapsed times shown in Fig. 6 were that of 
the SPH simulation and the times of the one-dimensional 
evolution were conveniently chosen to fit the 2D profiles. 

4.3 Implosion of a homogeneous capsule 

Probably the most unfavourable case to numerical simula- 
tion is that of a strong implosion as it currently appears, 
for example, in IGF studies. Thus, our third test deals with 
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Table 1. Main features of test models described in Sections 4.2 and 4.3. Conservation of momentum is given 
by the center of mass displacement Af = [{f (4) — fo) — S^tj and AS = |^(2(t) — zq) — divided by the 

radius of the configuration at that time (columns 5 and 6) . Momentum and energy conservation correspond 
to the last calculated model shown in Figs. 6, 7, 8 and 9. 
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the implosion of a homogeneous spherical capsule of size 
Rq = 1 cm and density po = 1 g-cm-'', induced by the 
(artificial) ablation of its surface. The ablation of the cap- 
sule was triggered by the instantaneous deposition of en- 
ergy in the outermost layers of the capsule, which had their 
internal energy increased in the same proportion. The en- 
ergy deposition profile was taken linear from s = 0.8 cm to 
s = R^) = 1 cm so that the internal energy at _Ro = 1 cm 
was a factor 10* higher than that at s=0.8 cm. Below s=0.8 
cm a flat profile of the internal energy was assumed. The 
rocket effect caused by the evaporation of the surface lay- 
ers forms a strong shock wave which compresses the interior 
of the capsule. The convergence of the shock at the center 
of the sphere increases the temperature and the density in 
a large factor, much larger than that obtained in the pre- 
vious section dealing with the collapse of a solar-like poly- 
trope. The main features of the model are shown in Table 
1 and Fig. 7. The profiles of density and radial velocity at 
several times are depicted in Fig. 7 which also shows the 
profiles obtained using a ID-lagrangian hydrocode (contin- 



uum lines) with the same physics and identical initial con- 
ditions. The evolution of the capsule can be summarized as 
follows. Soon after the initial energy deposition a pair of 
shock waves moving in opposite directions appear (profiles 
at t=0.0044 s in Fig. 7). As the reverse shock approaches 
the center it becomes stronger owing to the spherical con- 
vergence (profiles at t=0.0087 s in Fig. 7). Once the max- 
imum compression point is reached, pmax = 31 g.cm-^at 
t= 0.0111 s, the wave reflects. When t=0.0150 s the density 
peak has already dropped to 7 g.cm-^and most of the mate- 
rial of the capsule is expanding homologously. At t=0.0264 
s the reflected wave reaches the initial radius of the cap- 
sule, Rq = 1 cm. At that time the material of the capsule is 
rather diluted p(r) « 1 g.cm-^and the radial velocity pro- 
file consist of two homologously expanding zones separated 
by a transition region at s ~ 1 cm. At the last calculated 
time, t= 0.0264 s, the outermost layer of the sphere has 
expanded until s = 10 cm, ten times the original size of the 
capsule. As we can see the 2D-simulation match quite well 
the ID results. According to Table 1 the central density at 
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Figure 5. Density, gravity and pressure gradient profiles of the polytrope. Upper-left and right: density 
profile without and with axis corrections respectively. Bottom-left and right: the same but for the pressure 
gradient. The absolute value of gravity computed using the gradient of the gravitational potential calculated 
using equation (45) is marked with full triangles (bottom-right). All mass points of the polytrope have been 
represented. Neglecting axis corrections leads to a much larger dispersion in the profiles. 



the moment of maximum compression is almost the same 
in both calculations. The conservation of momentum is very 
good, close to machine precision (columns 5, 6 of Table 1), 
while numerical energy losses remained lesser than 1% (last 
column in Table 1). However, an inspection of Fig. 7 reveals 
that the spherical symmetry was not so well preserved as 
in the collapsing polytrope. Now the implosion of the cap- 
sule was more violent and the deceleration phase before the 
bounce was more intense, making easier the growth of in- 
stabilities. Even though the initial model had good spherical 
symmetry, the initial distribution of the particles in a regular 
lattice acts as a source of the so called hour-glass instability. 
Such instability is related to the existence of preferred di- 
rections along the grid through which the stress propagates. 
As in the previous test another point of conflict between the 
ID and the 2D calculations is that the 2D evolution is a bit 
delayed with respect its one-dimensional counterpart. For 
instance, the times at which the maximum central densities 



are achieved are t^o = 0.01095 s and t2D ~ 0.0111 s, so that 
the relative difference was around 1 — 2%. Such percent level 
of discrepancy remained approximately constant during the 
calculation. 

4.4 Wall heating shock: the Noh test 

The wall heating shock test, HHE l| 19871 ). was especially de- 
vised to analyze the performace of algorithms addressed to 
capture strong shocks. Basically the wall heating shock test 
consists in making a sphere or a cylinder implode by impos- 
ing a converging velocity field. For these geometries there 
is an analytical approach to the evolution of density and 
thermodynamical variables as a function of the initial con- 
ditions. The results of numerical codes can be compared with 
the analytical solution to seek the best method or to choose 
for optimal combination of parameters. It is well known that 
schemes which rely in artificial viscosity have difficulties to 
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Figure 6. Density profile (left) and radial velocity (right) during the implosion and further rebound of the polytrope at different times. 
The black continuum line is the profile calculated using a one-dimensional hydrocode of similar resolution. At t=1545 s the shock wave 
is breaking the surface of the star. In the figure all particles used during the calculation are shown. 




handle the wall heating shock test. The reason is that arti- 
ficial viscosity spreads the shock over a 3-4 computational 
cells, which induces an unphysical rise of internal energy 
ahead the shock. In the case of spherical or cylindrical ge- 
ometry the artificial increase in internal energy is magnified 
by the geometrical convergence of the shock. Therefore the 
wall shock problem represents a strong challenge for the ax- 
isymmetric SPH. Brookshaw (2003) carried out a similar 
test with the SPH code but far from the symmetry axis. In 
particular he modelled the impact of two supersonic streams 
of gas, obtaining good profiles for density and internal en- 



ergy except in a small region around the collision line. The 
density profile showed a dip at that region while a large 
spike was seen in the internal energy. Both, the dip and the 
spike are numerical artifacts which can be smoothed out by 
using an artificial heat conduction term to spread the excess 
of internal energy and reduce the error around the contact 
discontinuity. 

To check the performance of our code we have settled 
N=50334 particles in a square lattice. As in the previous 
tests the mass of the particles was conveniently crafted to 
reproduce an spherical homogeneous system with initial ra- 
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dius _ Ro — 1 cm. The initial conditions were taken as in 
iNohl j 19871 ): p(s,0) = 1 g.cm-^ s(s,0) = -1 cm.s'^ 
m(s, 0) = erg.g^^. The exact solution at time t=0.6 s 
for 7 = 5/3 is shown in Fig. 8 (dashed line). The analyt- 
ical profile is characterized by a constant post-shock state 
until distance s = 0.2 cm followed by a rapid decrease in 
density and internal energy. In the shocked zone density 
reach a constant value p = 64 g.cm~^while internal energy 
was u — 0.5 erg.g"^. The combination of such harsh initial 
conditions plus geometrical convergence leads to an strong 
implosion, even harder than that described in the ablated 
capsule test. In Fig. 8 (left) there are shown the spherically 
averaged profiles of density and internal energy obtained 
using SPH at t= 0.6 s. The error bars in the plot give the 
la dispersion of these variables with respect its mean value 
in the shell. As we can see the resulting profiles compare 
poorly with the analytical results in the shocked region. In 
addition the dispersion is high, especially at low radius, a 
clear signature for the presence of numerical noise. Close to 
the a:xis there is the typical artificial combination of a dip 
(in density) and spike (in internal energy). The maximum 
density value was p ~ 58 which is around 10% lower than 
the exact value. Such bad quantitative agreement was not 
unexpected because it is common to all hydrodynamic codes 
which use the artificial viscosity scheme. A way to improve 
the quality of the simulation is to allow for heat conduction 
in the hydrocode to remove the thermal energy spike and 
to sharpen the shock. Recipes to obtain an artificial con- 
ductivity coefficient in S PH to bett e r hand le th e wall shock 
proble m were given by iMonaghanI l|l992h and iBrookshawl 
l|2003l l. For this calculation we have adapted the recipe of 
Monaghan to the features of the axisymmetric SPH defining 
an artificial conductivity for the i-particle: 

ki = pij Cy^-hij (cij +4:\fj,ij\) (53) 

where Cv^^ is the symmetrized specific heat and |p,ij| is the 
artificial viscosity parameter given by equation (24). Accord- 
ing to equation (53) ki = kj, which can be used directly in 
equation (35) to compute the artificial heat flux. As shown 
in Fig. 8 (right) the inclusion of the artificial heat conduc- 
tion term leads to a significant improvement of the results. 
Not only the profound dip in density in the central region 
has been removed but the simulation shows much lesser dis- 
persion around the averaged values of density and internal 
energy. However the maximum peak in density still remains 
~ 10% below the theoretical value. In this respect the only 
way to improve the resu lts is sharpen i ng th e shock either by 
using adaptive kernels, lOwen et al.l (|l998h . ICabezon et al.l 
l|2008l . or by increasing the number of particles. 

4.5 Supersonic collision of two streams of gas 

Our last test is specifically addressed to check momen- 
tum conservation in a very anisotropic situation: the super- 
sonic collision along the Z-axis of two homogeneous spher- 
ical clouds of gas with different size and mass. The size 
and mass of the spheres are 7?i = 1 cm; Mi — 4n/3 g 
and R2 ~ S cm; A/2 ~ 36n g respectively so that their 
density is pi = p2 = 1 g.cm~^. The biggest sphere is at 
rest with its center located at (0, 0) cm while the smaller 
one is centered at (0, —5) cm moving upwards with velocity 



= 10 cm.s~^. The initial internal energy of both spheres 
was Ui = U2 — 10~^ erg.g"^ whereas EOS obeys an ideal gas 
law with 7 = 5/3. The corresponding initial Mach number 
was AI^ ~ 500 thus the impact is supersonic. The number of 
particles in each sphere was A''i — 5480 and = 50334 re- 
spectively. 

In Fig. 9 there are represented several snapshots show- 
ing the density evolution during the collision process. As can 
be seen the large mass contrast leads to the complete defor- 
mation of the smaller sphere which, in the end transfers 
most of its initial momentum to the larger bubble. Infor- 
mation about the evolution of the center of mass of each 
bubble as well as that of the whole system is provided in 
Fig. 10. According to Figs. 9 and 10 the collision history 
can be roughly divided in three stages: 1) For ^ f 55 0.5 s, 
the incoming smaller cloud deforms while a large fraction 
of its kinetic energy went into internal energy around the 
collision region. A shock wave was launched into the larger 
bubble, 2) Between 0.5 t 5C 1 s the total internal energy 
did not change so much. At t = 1 s the velocity of the center 
of mass of both structures was practically the same, 3) For 
t > 1 s the energy stored as internal energy is again restored 
to the system. At larger times the velocity of the smaller 
bubble became negative while the bigger cloud acquired a 
positive velocity to preserve total momentum. At our last 
time t = 1.4 s the interaction between both structures is 
coming to an end. Fig. 10 also shows that, in spite of the 
large changes in the velocity of each bubble, the velocity of 
the center of mass of the whole system remained unaltered. 

In the last row of Table 1 there are shown several magni- 
tudes related to momentum and energy conservation. Con- 
servation of linear momentum during the interaction was 
monitored through the displacement of the center of mass 
of the system. At the final time, t = 1.4 s, the deviation of 
the center of mass position with respect the value expected 
from v'^t was much larger than that shown in the other tests. 
Still the relative error in the position of the center of mass 
once the collision has practically ceased was ~ 10"^ An gu- 
lar momentum was very well preserved, almost to machine 
precision. Total energy was conserved to 10~^ relative to the 
initial kinetic energy of the system. 

In an attempt to understand the origin of the discrep- 
ancy relative to momentum conservation in the Z-direction 
we ran exactly the same model but this time symmetrizing 
equation (13) (multiplying the term Pjrj/r)^ by ff). There 
were no significant changes. We conclude that strict total 
momentum conservation in the Z-direction was not possible 
because of the interaction between real and reflected parti- 
cles. Such interaction took place in a small band around the 
symmetry axis, acting as an external force which modified 
momentum of real particles. However that force can not be 
balanced by an opposite force acting in the left semiplane 
because, in the Z-direction, reflected particles were obliged 
to move exactly in the same way real particles did. There- 
fore if strong directional anisotropics appear in the vertical 
displacements of the mass points some degree of violation 
in momentum conservation is unavoidable. Eventually mo- 
mentum conservation should improve as the number of par- 
ticles increase because the amount of mass settled in the 
axis vicinity is lower. 
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Figure 8. Results of the Noh test at t=0.6 s. Upper left and right: averaged density profiles without and 
with the artificial heat conduction. Bottom left and right: same but for internal energy. Dispersion around 
the averaged values are given by the error bars. 



5 CONCLUSIONS 



Despite the success of the smoothed particle hydrodynam- 
ics technique to handle gas-dynamics problems in three 
dimensions little effort has been invested to develop two- 
dimensional (axisymmetric) applications. This work intends 
to fill that gap by solving many of the problems often 
invoked in connection with SPH in cylindric coordinates. 
These are the treatment of the singularity axis, the han- 
dling of shock and thermal waves and, in many astrophysi- 
cal scenarios, an accurate procedure to calculate the gravita- 
tional force. Our general philosophy in developing the math- 
ematical formalism was to remain as close as possible to the 
cartesian scheme so that many of the results of the standard 



SPH in three dimensions can be extrapolated with minimum 
changes to the axisymmetric version. 

Starting from the fluid Euler equations given by 
iBrookshawl (j2003l ) we have obtained analytical corrections 
to those particles which move close to the singularity axis. 
These corrections appear as multiplicative factors to the dif- 
ferent terms of the fluid equations. Such multiplicative fac- 
tors become equal to one when r > 2h (Fig. 1) so that 
there is no need to calculate them beyond that distance 
and the method is computationally efficient. Once the ba- 
sic formalism was built we added several pieces of physics 
which make the scheme well suited to handle a large vari- 
ety of problems. First of all an extension of the 3D stan- 
dard artiflcial viscosity to the 2D-axisymmetric realm was 



Axisymmetric smoothed particle hydrodynamics with self-gravity 15 



"T 1 1 r 




J I I L_ 



J I L_ 



II .5 



10.5 



Figure 9. Density colour-map of the collision of two streams of gas at times t= s, t= 0.31 s, t= 0.69 s and t= 1.41 s. Axis units are 
in centimeters. 



devised, equations (23), (24) and (25). The artificial viscos- 
ity includes the convergence of the flow towards the sym- 
metry axis via linear and quadratic terms proportional to 
jtjj^j (^|-,gjjjg J. ^jjg distance to the symmetry axis). These 
terms arise because the diagonal part of the stress ten- 
sor in cylindric coordinates includes the divergence of the 
flow velocity. However the cartesian part of the divergence 
was already present in the standard formulation given by 
equations (22) and (23) giving rise to bulk and shear vis- 
cosity. Therefore only the axis converging part of velocity 
divergence, proportional to Vr/r has been added. The re- 
sulting viscous force H^'^ has two terms: n*^' (cartesian) 
and n'^-' (axis converging part) which are calculated using 
similar expressions, equations (23), (25) involving four con- 
stants ao,/3o (n'^') and ai,/3i (n'^-*). The axis converging 
part of viscosity will be of importance only for strong im- 
plosions. Although in the simulations reported in this paper 
we have taken ao = «i = 1 and /3o = /3i = 2 there could be 
other plausible combinations worth to explore. A similar ap- 
proach was used to built the conductive transport equation. 
Starting from the expression given bv .Brookshaw (,1985 ) we 
have added a new term that accounts for the divergence of 
the temperature gradient in the axis neighborhoods. Even 
though the resulting expression was not totally antisymmet- 
ric it led to a satisfactory energy conservation in the tests 
described in sections 4.1 and 4.3. Finally the set of equations 
was completed with the inclusion of self-gravity. In axisym- 
metric geometry it is better to rely in the direct, ring to ring, 
interaction to compute gravity. Such procedure, although in 



general computationally expensive, has several advantages: 
1) it does not rely in any sort of grid and gives the exact 
value of the gravitational force, 2) the calculations are feasi- 
ble for a moderate number of particles (i.e. around 5 10* in 
serial desktop machines) which often suffice in many 2D sim- 
ulations 3) the scheme is clear and extremely simple making 
it straightforward to parallelize. The algorithm devised to 
compute gravity was able to keep track the free-fall collapse 
of an homogeneous sphere as well as of an homogeneous ro- 
tating cylinder. Given an initial amount of rotation the evo- 
lution of any point of the cylinder was handled by imposing 
angular momentum conservation around the symmetry axis. 

Five test cases aimed at validating the computational 
scheme were considered. Each of them intended to check a 
particular physics item: heat conduction (Sects. 4.1 and 4.4), 
shock waves (Sect. 4.3 and 4.4) and gravity (Sects. 3.3.1 and 
4.2). Momentum conservation was specially checked in Sect. 
4.5. Because for the most part the initial configurations had 
spherical symmetry the preservation of that symmetry dur- 
ing the collapse and further expansion of the system was 
taken as an indicator of the general ability of the scheme to 
handle large changes in the spatial scale. It was also possi- 
ble to make detailed comparisons of the results with those 
obtained either from analytical calculations, as in the case 
of the thermal wave propagation, or from one-dimensional 
simulations carried out using standard lagrangian hydrody- 
namics. The results, obtained using several dozen thousand 
particles, were in good agreement with the analytical expec- 
tations and to the ID simulations. Additionally there was 
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Figure 10. Evolution of Vz of the center of mass of eacli bubble 
and of the system. Subscripts 1 and 2 refers to the smaller and 
larger bubble respectively. Evolution of total kinetic and internal 
energy is also given. 



an almost perfect momentum conservation while numerical 
energy losses always remained lesser than 1%. The spher- 
ical symmetry was well preserved although a slight devia- 
tion was observed in the capsule implosiou test, partially 
due to the so called hour-glass instability (owing to the ini- 
tial setting of the particles in a lattice) and, despite the 
corrective terms, to the effect introduced by the singularity 
axis onto the approaching particles. A similar but probably 
more demanding test was the wall heating shock problem, 
dealing to the spherically symmetric implosion of a super- 
sonic stream of gas. Although the results are not as good 
as those obtained with other methods which do not use ar- 
tificial viscosity, they are comparable with one-dimensional 
calculations with similar resolution relying in artificial vis- 
cosity schemes. However the results improved when artificial 
heat conduction, calculated using equations (35) and (53), 
was switched-on. 

As discussed in Sect. 4.5 a weak point of the formu- 
lation is that the inclusion of refiective particles in the 
scheme could degrade total momentum conservation in the 
Z-direction. However good momentum conservation, much 
better than energy conservation for instance, is expected 
in those systems in which the velocity field is not too 
anisotropic. If strong momentum exchange along the Z-axis 
is expected, as in the test described in Sect. 4.5, then mo- 
mentum will be preserved to a similar level as total energy. 

As a conclusion we can say that the formulation of the 
axisymmetric SPH given in this paper is a solid tool to carry 
out simulations using that particle method. However there 
are still a number of unresolved issues which deserve fur- 
ther development. One of them is how to build good ini- 
tial models without settling particles of different mass in 
regular lattices. This is important because ordered lattices 
introduce spurious instabilities in the system and the mix- 
ing of particles with very different masses could lead to nu- 



merical artifacts. Another point of difficulty has to do with 
artificial viscosity, because it introduces too much shear vis- 
cosity in the system damping the natural development of 
hydrodynamical instabilities. The new term H^^' given by 
equation (25) of artificial viscosity comes from the diagonal 
of the stress tensor so its contribution to shear viscosity is 
probably lower than that of H^^^ . In amy case the artificial 
viscosity formulation given in this work is so close to the 
standard formulation that it could benefit from future ad- 
vances in the much better studied three-dimensional SPH. 
Finally, the simple test about the implosion of a rotating 
cylinder indicates that axisymmetric SPH is able to handle 
rotating structures but more work needs to be done to in- 
corporate angular momentum transport into the numerical 
scheme. 
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5.1 Appendices 

APPENDIX A: CORRECTION FACTORS 
CLOSE TO Z-AXIS 

First, wc (Icrivc the factor /i affecting the 2D-dcnsity q. 
Close to the Z-axis symmetry enforces the 3D-density p to 
have a maximum or a minimum. Thus we can safely assume 
p — po in the axis vicinity. If wc make the one-dimensional 
SPH estimation of the averaged density along the r-axis: 



)dr' 



(Al) 



where W^^ is the cubic-spline kernel in ID: 



^ 



l<Ur^2 (A2) 



Ur > 2 



being Ur = \{r — r')\/h. If we take ri{r') = 2-n\r'\po, which is 
only valid for r' — » 0: 



(r7Wr)i = ^^mjVrjW{\si - Sj|,/ii) x /j 



Similarly the correction to be applied to riVz is: 



< V^z >= 



/i(C) 



and its corresponding discrete expression is: 



{'nvz)i = I rrijVz 



(A8) 



(A9) 



(AlO) 



A plot of the factors /i(C) and /2(C) and their radial 
derivatives is given in Fig. 1. The factor /i(C) is exactly one 
for r = 2h and goes to zero when r/h 0, as expected. 
However, it can be seen that factor /2(??) is close but not 
exactly one at r = 2h and its slope is not totally flat. This is 
because the integral expression in equation (A6) is quadratic 
in the radial coordinate r thus the interpolation does not 
give the exact value of < r]Vr >■ 



<7){r) >=2%po / l'r'lH/'"( 



h 



^-)dr' 



(A3) 



Using the cubic-spline the right side of equation (A3) can 
be integrated to give: 



< r){r) >= 



n 



MO 



(A4) 



where C,= r/h and rj = 2nrp is the corrected density (here- 
after we put a hat over any corrected quantity). The cor- 
rection factor /i(C) is given by equation (5). The density 
in brackets is what SPH computes using summatories in- 
stead of integrals. Thus, adding the z-coordinate and using 
W^'^ the corrected density can be evaluated using: 

?(r) =< 7? > /i(C) = |^f^m^I¥^°(|ri - rjU^j x /^AS) 

A similar procedure can be used to get the adequate ex- 
pressions for riVr and rjVz, needed to compute the temporal 
evolution of density, equation (15). Symmetry enforces the 
radial velocity to vanish at the symmetry axis, thus close to 
that axis Vr = Cr and: 



00 

<nvr >= piz)2nK^°C j \r'\r'Wr{ur)dr' 



(A6) 



Again the integral on the right admits an analytical 
solution for the cubic-spline kernel. After a little algebra the 
following expression is obtained: 



< riVr >= 



MO 



(A7) 



where /2(C) is given by equation (16). The expression inside 
the brackets is what the numerical code calculates using 
summatories instead of integrals. Thus, the corrected rjur, 
value of that magnitude is: 



